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We study the global bending modes of a thin annular disc subject to both an internally 
generated magnetic field and a magnetic field due to a dipole embedded in the central 
star with axis aligned with the disc rotation axis. When there is a significant inner 
region of the disc corotating with the star, we find spectra of unstable bending modes. 
These may lead to elevation of the disc above the original symmetry plane facilitating 
accretion along the magnetospheric field lines. The resulting non-axisymmetric disc 
configuration may result in the creation of hot spots on the stellar surface and the 
periodic photometric variations observed in many classical T Tauri stars (CTTS). 
Time-dependent behaviour may occur including the shadowing of the central source 
in magnetic accretors even when the dipole and rotation axes are aligned. 
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1 INTRODUCTION 

Situations in which a thin accretion disc is threaded by a strong poloidal magnetic field are of interest in different astrophysical 
contexts relating to accreting objects such as neutron stars or young stellar objects. This situation may result when the disc 
interacts with a magnetic field produced by a dipole embedded in the central star such that the dipole field lines penetrate 
the disc (Ghosh & Lamb 1978, Ghosh & Lamb 1979, Campbell 1987, Camenzind 1990, Konigl 1991). 

In addition open poloidal field lines may be advected inwards by the accreting matter (Lubow, Papaloizou & Pringie 1994, 
Reyez-Ruiz & Stepinski 1996, Agapitou & Papaloizou 1996) and be associated with a centrifugally driven wind (see Konigl 1993 
for a review). The observed correlation between mass accretion rate and mass outflow rate in T Tauri stars (TTS) supports 
the idea of the accretion disc as the underlying source of the outflows emanating from young stars (Cabrit et al. 1990). 
The short-term photometric variability evident in many CTTS has been attributed to the presence of both dark and hot 
spots that cover a part of the stellar surface (Bouvier 1994, Bouvier et al. 1995). The latter have been explained as arising 
from shocks formed close to the stellar surface resulting from non-axisymmetric accretion along stellar magnetic field lines. 
Such magnetospheric accretion is also invoked to explain the large infall velocities inferred from emission lines (Calvet & 
Hartmann 1992, Edwards et al. 1994, Hartmann, Hewett & Calvet 1994), the infrared colours of TTS (Kenyon, Yi & Hart- 
mann 1996) and the outbursts of EX Lupi (Lehmann, Reipurth & Brandner 1995). Most of the variability of TTS cannot be 
accounted for unless there is some time-dependence in the magnetic fleld and(or) the accretion flow (Bouvier et al. 1995). 
The presence of an accretion disc around many TTS seems also to be linked to their low rotational velocities and the kind of 
activity observed on the stellar surface (Edwards et al. 1993, Montmerle et al. 1993). 

Under some conditions, the interaction of a strong stellar dipole fleld with a surrounding low mass disc prevents accretion. 
This is because the magnetic stresses exerted on the disc, external to the radius where the star corotates with the local disc 
material which has near-Keplerian rotation, act so as to transfer angular momentum to the disc and spin down the star. The 
accretion disc is then truncated at a radius where the viscous and magnetic torques balance. Models of this type have been 
developed to account for the low rotational velocities of TTS provided that the accretion disc is not dissipated too early in 
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the star's lifetime (Cameron & Campbell 1993, Yi 1994, Ghosh 1995, Armitage & Clarke 1996). We note though that non- 
magnetized stars could also experience spin-down while they accrete mass. Paczynski (1991) and Popham & Narayan (1991) 
have shown that when a star rotates close to its breakup speed the accreted specific angular momentum decreases and it even 
attains negative values for large enough values of the stellar rotation speed. 

The situation of no or low mass accretion together with stellar spin down occurs when the field is strong and the disc mass is 
low giving rise to a small viscous angular momentum flux. In this paper we shall consider the converse situation when the disc 
is massive enough to penetrate through to the corotation region so that accretion becomes possible. As in previous studies, we 
assume the existence of a dipole field embedded in the central star. The present observational evidence cannot rule out such 
a coherent field structure (Montmerle et al. 1994) . There is evidence from numerical calculation of non-linear stellar dynamos 
that a steady dipole mode is the most easily excited one (Brandenburg, Tuominen & Moss 1989). The possibility of a fossil 
field has also been proposed (Tayler 1987). 

In this paper we adopt the model of Spruit & Taam (1990) which is such that when the disc is able to reach interior to the 
corotation radius, the inner parts which contain field lines connected to the central star, corotate with it. In that case, simple 
thin disc models in which the magnetic field plays an important part in supporting the inner corotating region against gravity 
can be constructed when the axis of the dipole and the disc angular momentum axis are aligned (Spruit & Taam 1990, 1993). 
Considering the inner corotating disc to be part of the magnetosphere, these models have the fastness parameter, or ratio 
of stellar to inner differentially rotating disc angular velocity, close to unity. Magnetic support may also be important in the 
outer, differentially rotating part of the disc if an inwardly advected field becomes strong, as would be the case in the presence 
of a strong wind removing most of the disc angular momentum (see for example Konigl 1989). 

The importance of the issue of the stability of these disc models has been stressed by Spruit & Taam (1990) who considered 
the role of interchange instabilities in enabling matter to migrate inwards in the inner corotating region of the disc until direct 
particle motion confined to equilibrium vacuum field lines becomes possible. However, note that such particle motion along 
vacuum field lines may not be representative of the plasma fiow that may occur along field lines because the plasma may in 
principle contain significant currents that disturb the original vacuum field lines. A study of the possible flow along field lines 
accordingly requires a full MHD treatment. The local stability of the outer differentially rotating disc to interchange modes 
has been considered by Spruit, Stehle & Papaloizou (1995). 

In this paper, we study the global stability of a thin magnetized accretion disc to both axisymmetric and non-axisymmetric 
disturbances perpendicular to its plane (bending modes). These are the thin disc limit of modes with density perturbation 
having odd symmetry with respect to reflection in the mid-plane, in contrast to the interchange modes which have even 
symmetry. At equilibrium the disc is permeated by both an internally produced poloidal magnetic fleld and an external dipole 
field. We here limit consideration to the case when the vertical component of the field in the inner disc does not change sign. 
Bending modes are of potential interest because an instability may lead to elevation of the disc mid-plane above the original 
symmetry plane leading to facilitated motion along field lines connected to the star as well as time-dependent shadowing 
of the central source. In the context of neutron stars this was also pointed out by Spruit & Taam (1990). The accretion 
along magnetic field lines derived from a non-axisymmetric disc would result in the production of hot spots on the stellar 
surface and modulation of the power output. This could lead to a time-dependent accretion flow even in the aligned dipole 
case and it may account for the irregular variability of CTTS without the need to invoke a variable magnetic fleld. Finally, 
non-axisymmetric modes with azimuthal mode number m = 1 are related to disc precession (Papaloizou & Terquem 1995) 
and are of potential interest with regard to precessing jets. 

In section |2| we describe the thin equilibrium disc models. In section ^ we give the perturbation equations for linear modes 
under the ab initio assumption of a razor-thin disc. We derive the local dispersion relation which indicates instability in the 
inner corotating parts of the disc if they are extensive enough. We further derive variational principles for the axisymmetric 
modes showing that the unstable modes do occur in discs with flnite albeit small thickness. We give a simple physical picture 
of the instability showing how it originates as an unstable interaction between the central dipole and current loops in the 
disc. In section ^ we describe the speciflc equilibrium disc models that we consider and our numerical calculations. They give 
spectra of axisymmetric and non-axisymmetric unstable modes conflned to the disc inner regions. These results are presented 
in section ^. Finally in section ^ we discuss possible consequences of our results. 



2 EQUILIBRIUM DISCS 

We consider thin disc conflgurations with an axisymmetric poloidal magnetic field B — {Br,0, B^). Here we use cylindrical 
polar coordinates {r,(f),z). The field is described by a fiux function i/) such that 

-Ihijj 1 dip , . 

Br ^ — ^andB^^- — . (1) 

r oz r or 

For this field the current density j — (0,^,^,0). For an infinitesimally thin disc it is convenient to work with the vertically 
integrated azimuthal component of the current density J, where 
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By integrating the azimuthal component of Ampere's law through the disc, we obtain 

B+ = J/2, (2) 

where -B^ denotes the radial component of the magnetic field on the upper surface of the disc. Here we have set the magnetic 
permeability, fio, to unity. We can recover the equations in MKSA units by replacing B by B/y^po and J by J/^/JTq. Br is 
antisymmetric with respect to reflection in the disc mid-plane so that its value on the lower surface of the disc is B^ ~ —B^. 
Thus Br changes significantly on passing through the disc in contrast to Bz which, as implied by the condition V ■ S = 0, 
changes negligibly on passing through an inflnitesimally thin disc. 



2.1 Force balance 

The vertical integration of the radial component of the momentum equation yields the condition for radial equilibrium as 

E— ^Ern^ + JB,, (3) 
or 

where E is the surface density, fl is the disc angular velocity, and $ is the gravitational potential. 

The above formalism neglects the radial pressure force. An inflnitesimally thin disc can in practice be considered to be 
signiflcantly thinner than Cs/Q (cs being the sound speed). The vertical equilibrium then implies that the mid-plane pressure 
P ~ {Br) 1 since the magnetic squeezing of the disc overwhelms its tidal confinement. Thin disc equilibrium configurations 
of strongly magnetized discs with the above properties, which take account of the vertical structure, have been constructed by 
Ogilvie (1997). In order to neglect pressure forces in the radial direction, we require, assuming B^ ~ Bz, that H <^ Lr, where 
H is the semi-thickness of the disc, and Lr is the scale length of variation in the radial direction. Thus the approximation 
scheme becomes better for thinner discs. We consider equilibria for which the gravity is due to a central point mass, M, such 
that 

GM 



Equilibrium models may be constructed with the surface density, E, and integrated current density, J, being specified as 
arbitrary functions of r. The flux function with no external sources is then given by (see Lubow et al. 1994) 



^(r,.) = -/ / J(rOcos(^')/d^'dr' 

47t /o /„ y/r^ + r'2 - 2rr' cos(v9') + z'^ 



where the disc is presumed to have inner and outer boundary radii Ri and Ro respectively, the latter possibly being infinite. 
To the above internally produced flux, we may add a contribution due to external sources, i/'ext- When the external source is 
a dipole at the origin, 

^- = -sr(i?i)-^^, (5) 

(r-^ + z^) ' 

where Bz^*{Ri) is the external vertical field at the disc inner boundary. 

In the presence of a central dipole, some of the field lines which cross the inner regions of the disc may join to the dipole in 
the centre. Further out, field lines may be open in the case of an infinite disc or close before the outer boundary when the 
disc is finite. Field lines with these properties appropriate to an equilibrium configuration are illustrated in Fig. |^ below. The 
sign of the azimuthal current density in this and all other configurations we consider here is such that the sign of the vertical 
field in the inner regions of the disc does not change. That is there is no X point. This is the situation naturally expected if 
dipole field lines diffuse into the disc. 

For physical consistency, field lines joining the central dipole should be in a state of isorotation at constant angular velocity. 



so E should be specified accordingly (see section O). The condition of isorotation means that the magnetic field must make 
an increasingly important contribution to support the fluid against gravity as r decreases. Spruit & Taam (1990) have argued 
that material moving inwards from the outer disc due to angular momentum transport processes occurring in accretion discs 
(see Papaloizou & Lin 1995 for a review) can migrate inwards into the isorotating region due to the action of interchange 
instabilities. This may produce a magnetically dominated isorotating thin disc if the material remains cool. In fact in order 
to establish an inner corotating part of the disc, magnetic support against gravity should not be large at the outer corotation 
radius. Assuming 5+ ~ Bz, this requirement gives at that radius 

Bl « 
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Using M = 2nrvrT., Vr being the radial velocity, and the viscous inflow rate Vr — —vjr, v being the kinematic viscosity, we 
obtain 

Bl « 



Thus, as indicated above, for flxed stellar properties and disc viscosity, establishment of an inner corotating region is favoured 
at large accretion rates, \M\, and accordingly large disc masses. 

Open field lines in the outer disc may in principle rotate at any angular velocity. Field lines that close in the outer regions of 
a finite disc may be opened if there are additional external currents which could be produced by, for example, a wind. 



3 PERTURBATION EQUATIONS 

We consider linear perturbations of the equilibrium configurations with a Lagrangian displacement { which, for a razor-thin 
disc, has the form 

i = (0,0,^.). 

The only non-negligible component is the vertical one which is independent of z. This displacement belongs to a class such that 
is even, while the other components are odd with respect to refiection in the disc mid-plane. These are thus bending modes. 
We may also assume that the (^-dependence of the perturbations is through a factor exp(im(/3), m denoting the azimuthal 
mode number. From now on this factor will be taken as read and will be dropped from the perturbations. The Eulerian 
perturbations of the various quantities are denoted by a prime. 

The perturbation of the magnetic field interior to the disc, B' , is related to ^ by the integration of the linearized induction 
equation with respect to time 

B' = {Bl,B'^,B'z) = Vx(ixB) (6) 
The non-zero components of B' take the form 

B; = -e.^,andB-i^(i|ii^ (7) 

az r or 

where B'^ is antisymmetric with respect to refiection in the disc mid-plane and B'^ is symmetric. 

The vertical component of the perturbed field in the disc must be matched to the vertical component of the perturbed vacuum 
field exterior to the disc. The perturbed vacuum field may be taken to be a potential field. This is the case even when the disc 
takes the form of an annulus making the vacuum multiply connected, because the symmetry properties of the magnetic field 
perturbation make it circulation-free. On the upper disc surface we therefore have 9$!y[/5« = B'^'^ , where is the value 
of the vertical field perturbation just outside the disc surface and $m is the magnetic potential associated with the external 
field perturbation. Continuity of the vertical field component at the upper disc surface implies that $5^ can be found using 
(^). Thus we obtain that on the upper surface 

1 9(ri?+g.) 

^ — ^ s • y°) 

oz r or 

The corresponding equation with + ^ — applies on the lower surface. Finding $5^j is entirely analogous to finding the 
gravitational potential ^I" due to a disc surface density distribution (see Tagger et al. 1990, Spruit et al. 1995). If $m is taken 
to be equivalent to ^t, the appropriate surface density is equivalent to B'^ /{2-kG), where G is the gravitational constant. Thus 
$M may be written in the form of a Poisson integral 

^'^ = _i r B'+{r')cos{m^'yAr'A^' 



JiJi Jo \Jr''^ + r^ - 2rr' cos(v3') + z^ 

The radial component of the magnetic field perturbation on either the upper or lower surfaces of the disc, just outside the 
disc, is then given by 

B'+ = B'- = f^) . (10) 

3.1 Vertical component of the equation of motion 

In general, the vertical component of the Lorentz force per unit volume is 

F.^-i'-f+B.'-^. (11) 
2 oz or 
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Perturbing and integrating this vertically through the disc gives 

F'.dz^ ~2B+B'+ -2f,B+^ (12) 

or 

where we have assumed that is almost independent of z in the disc. The perturbed vertically integrated a-component of 
the equation of motion is 

where D denotes the convective derivative. 

The coefficients of equation (^|) are independent of t. We can therefore look for solutions in the form of normal modes. In this 
case the time-dependence of the perturbed quantities is taken to be through a factor exp{iat), where a is the eigenfrequency 
of the mode. Using this together with the fact that for a point mass potential 

where Qk is the keplerian angular velocity, equation (|l3[) becomes the normal mode equation 



((J + mD,) — $1k 



2B+ dB, 
E ~97 



«.-^f5^j__^ (») 



Equation (^) together with and (Q) constitutes a linear eigenvalue problem with a as the eigenvalue and as the 
eigenfunction. 



3.2 Local Dispersion Relation 

We can derive a local dispersion relation from ( p^ ) by adopting perturbations of the form cx exp(ifer), where k is the radial 
wavenumber, assumed to be ^ |m|/r. We comment that because the evaluation of $m is equivalent to calculation of the 
gravitational potential due to a surface density B'^ /{2nG), the situation here is closely analogous to that for bending modes 
in a self-gravitating disk (see Hunter & Toomre 1969, Shu 1984). The integral in (^) can thus be calculated using the WKB 
approximation to give 

4>M = -B'J\k\ = -ikB+^z/\k\ 

where $m is now the amplitude of the mode with radial wavenumber k. The local dispersion relation derived from (^^ is 
then 

{(j + mQ.) =Q.yi + —^— h ' ' fc . (15) 

2j or L 

Instability ensues on the existence of at least one mode with growth rate \{a-\-mQ,) > 0, for which we require ((T + mfJ)^ < 0. As 
the last term on the right-hand side of ( p^ ) is positive definite and therefore stabilising, the condition for instability becomes 

2 , 25+ dB^ 

f^K + — <0. (16) 
L or 



When the gravity is due to a central point mass, (^ gives 



^(J7|-f7^) = ^^, (17) 



so that (|l|) may be expressed in the equivalent form 



+ ^^<0. (18) 



(nl - Q?) Bz 9r 

We comment that in the non-rotating case (Q. = 0), ( p^ becomes the same condition as that given by Wu (1987) and 
Lepeltier & Aly (1996) who considered non-rotating current sheets. The condition ( p^ becomes that given by Anzer (1969) 
and Spruit & Taam (1990) provided one sets SIk = 0. This is because these authors omitted gravitational restoring forces in 
the direction perpendicular to the current sheet. 

The local criterion for stability becomes satisfied when the disc is being primarily supported by an external field. For an 
external dipole 

r QBz ^ 
Bz 9r 
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Then (|l^) will be satisfied when the magnetic field provides enough support against gravity so that > 30^/2. The latter is 
satisfied in the interior regions of the disc where the field lines link to the central dipole and Q is constant with 57k increasing 
inwards. 

We note that the condition above has been obtained for small |fcj , which is out of the domain of validity of a local approximation. 
However, the term in (jl^) involving |fc| is proportional to (-B^)^ and it therefore becomes of decreasing importance for a 
uniformly rotating magnetically supported region as this extends inwards towards regions of large B^, with the result that 
instability must eventually ensue. However, the precise details of onset require explicit calculation. 



3.3 Axisymmetric Modes 

Rigorous global criteria may be obtained in the case of axisymmetric modes through the existence of variational principles. 
In this case, for m — 0, ( p^ ) can be written in the form 

a^?. = (19) 
where the operator O is self-adjoint in that for arbitrary eigenfunctions and -q^ the following equality is satisfied 

In this case, a sufficient condition for instability is that, for any 



Ri J R 



1 



TT 



Ra fRo r^TT 



R: J Ri Jo 



{B'+{r)y B'+{r') cos{mip')rr'drdr'd^' 



+ - / / 'r-^ , . < 0- (20) 



^ j,/2 _|_ J.2 _ 2rr' cos(iy9') 



On insertion of suitable local trial functions, this gives the same condition as (|l^) . 
3.4 Relation to thick disc analysis 

We here note that the above variational principle may also be derived from the general variational principle of Papaloizou & 
Szuszkiewicz (1992) for stability to adiabatic perturbations of a general difi'erentially rotating equilibrium with a purely 
poloidal magnetic field, when the thin disc limit is taken. 

This establishes that the results are not an artefact of the use of the razor-thin disc approximation and vertically averaged 
equations from the outset. A sufficient condition for stability to axisymmetric modes is that for any trial 

C-L{i)dV<0 (21) 

where i is a linear operator, defined in Papaloizou & Szuszkiewicz (1992) and below, similar in principle to O but which 
depends on both r and z. The integral is now a volume integral (see below). 

The condition (|l]) is also necessary for stability to density perturbations with odd symmetry with respect to reflection in the 
equatorial plane such that {£,r,^^,^z) {—^r,—^ip,^z)- This corresponds to modes of the symmetry type considered in this 
paper. For these we have instability if 

- y r ■ = J (^-^ + pQii, r ) - ^ ■ v(vy)* + r ■ v^') + dv < o. (22) 

Here, p, Fi and ^' denote the density, specific heat ratio and magnetic flux perturbation respectively. For a fluid with smoothly 
vanishing density and pressure at the boundary all integrals, other than the last, are taken over the fluid volume V. The 
last integral which is related to the magnetic energy associated with the perturbation must be taken over the whole of space 
excluding any perfect conductors. The pressure perturbation is given by 

P' = -FiPV-^-C- VP. 

The quadratic form Q(^,^*) is given by 

Qite) = r{e ■ vr)(i . v^f) - r ■ ^1 ■ (J^ ~^)~ ■ ^ {^) ■ (^s) 

In our case we adopt ^ — (0, 0, ^z) with being real and depending only on r. For this trial function, use of tp' — — ^ ■ Vi/i 
inside the ideal fluid, together with vertical hydrostatic equilibrium, gives 
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Here all integrals except the last are taken over the fluid disc while the last is taken over the vacuum outside. 
We see that there is a close correspondence between (|^) and which applies to the razor-thin disc. Each of the integrals 
in (^) apart from the second, which becomes vanishingly small compared to the first in the thin disc limit, approaches the 
corresponding integral in (ISB). Thus the instabilities we find are not due to the ah initio assumption of a razor-thin disc. 



3.5 A simple picture of the instability 

Ifere we show how the condition for instability ( [l6| ) can be found from a simple argument relating to the interaction of a 
current loop with the central dipole. For the dipole fleld interior to the loop to be of the same orientation as the field produced 
by the loop itself (no X-point), the dipole moment must have the opposite sign to the azimuthal current density which 
produces an unstable interaction. A current loop at (r, z) with total current I produces a magnetic fleld at the origin with a 
vertical component 

^0 - 2(,2 f ^2)3/2- (25) 

The energy of an anti-aligned dipole of dipole moment ^ situated at the centre of the loop is 

^ = ''^°= 2(^2+^2)3/2- 

If we suppose that the loop lies initially at 2 = 0, and that is then displaced vertically up to z, where z is small compared to 
r, the change on the energy is 

4j.3 

To obtain the total change in energy we must add the change in gravitational potential energy msS7|;2^/2, where ms is the 
loop mass. Writing ms = 27tErdr, and / = Jdr, we obtain for the total energy change 

5W = + ^k] TtErz^dr. 

Using (^) and the expression for the external vertical fleld produced by the dipole, 3°]^^ — fi/{4nr^), we obtain 

5W=( + n^rz'dr. 
\ E or J 

The natural condition for instability is that energy be released on making the displacement, or in our case SW < 0. This gives 
the same condition as (|l^) if we replace Bz with the external vertical field which is the one that dominates under conditions 
of strong magnetic support. 

Thus we see that the generic instability can be understood in terms of an unstable interaction between the disc current and 
the central dipole. In the simple example described above, when the system passes through marginal stability because of say 
an increasing dipole moment, a bifurcation occurs such that there is a new stable equilibrium for the current loop lying off 
the original mid-plane. This suggests that in the case of a continuous disc, a warped structure may be taken-up in which the 
inner regions are elevated above or below the symmetry plane, facilitating eventual motion towards the central object along 
field lines. We now go on to discuss numerical calculations of normal modes for some specific models. 



4 NUMERICAL CALCULATIONS 

4.1 Equilibrium disc current and flux function 

In order to solve (|l^) we need to specify functional forms for E and J. The radial and vertical components of the magnetic 
field are then determined from ip which is calculated using (^. To evaluate the integral in the right hand side of (^ we use 
the method outlined in Lubow et al. (1994). The region between Ri and Ro is divided into Nr equally spaced grid points 
with separation dr. The surface current density J is discretised in the form of concentric ring currents I{ri) — J{ri)dr where 
subscripts i (and j later) denote values calculated at the i*** (j*^) grid point. Equation (^) is then approximated by 

V'(n) = ^n^7^,,J(r,)r, (26) 
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where Kij is the discretised form of the integral with respect to ip in (^); this can be written in terms of elliptic integrals (see 
Jackson 1975). We thus write 

4 [(2-fc')£i(fc")-2S2(fe')] 



y^n^ +rj'2 + 2ri rj k 
where Ei{k^) and E2{k^) are elliptic integrals of the first and second kind respectively, and 

In the numerical calculations performed here, the surface current density J is, to within a scaling factor, taken to be given by 

J = fM)~ fj{R, + Ar) 

where 

/.T(n) = c^V^ [exp(c2 ■ nc{u ~ n^id)V^o') + {rJR,y^-"T^'^\ (27) 

Both J and E can be taken to vanish at some point simultaneously in the same manner so as to retain a finite Lorentz force. 
We have chosen this point to be at a fictitious additional grid point at i? = i?o + dr for numerical convenience. 
The values of the constants Ci . . . C3 , r^i^ and nc are : 

Cl =10 C2 = 0.1 C3 = 2 

nc — A rniid ~ 0.9-Ro 

The radial distribution of 1/) is calculated from equation (|6|). The corresponding field lines in the z"*" quarter of the disc are 
the contours of i/)(r, z) = const. To the above internally produced magnetic flux we add a contribution from a dipole at the 
origin given by (^. Since the disc is thin, we shall neglect the radial magnetic fleld produced by this dipole and consider only 
the vertical component. 

4.2 Dimensionless variables 

We define a modified epicyclic frequency Km such that 

-Sn=fi^ + ^^ (28) 

L or 

and consider the dimensionless variables x — r/Ro, R = Km 

E = E/So, with 0,(1 — flK{Ro) and Bq and Eq being some fiducial values of the magnetic field and surface mass density 
respectively. The normal mode equation (|l^) can then be written in the dimensionless form 

(K^-(a + mn)^)c. = -/3oH^^ (29) 

were = $5y[/(_Bo^?o) and /3o is a constant which measures the relative strength of the magnetic and centrifugal forces (the 
larger /3o, the larger the magnetic support): 

BqRo 

The dimensionless external magnetic flux in the disc {z — 0) takes the form -i/jcxt ~ —"tpo/^- Equilibrium models which have 
scaled specifled proflles for disc current and surface density are parametrised by the two parameters, tf>o which measures the 
ratio of dipole to disc flelds and /3o which can be thought of scaling the disc surface density to provide a desired degree of 
magnetic support. 

In section ^ we present the results of normal mode calculations using equilibrium models with two different values of ipo, 
models of type I with t/jq = 0.03 and models of type II with tpo — 0.06. Various values of /3o have been considered. We plot the 
dimensionless components of the magnetic fleld produced by the disc currents in Fig. |l| The magnitudes of the dipole fleld 
at the innermost radius are given for models of type I and II for comparison. The contours deflned by i/)ext + tp — const are 
plotted in Fig. |^(a) for models of type I and in Fig. ^b) for models of type II. 

4.3 Surface density and Angular velocity 

The surface density and angular velocity are calculated in the following way. The magnetic support, Sm, is deflned as the ratio 
of magnetic to centrifugal forces in the condition for radial equilibrium: 
, , 2|B+BJa;2 
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Figure 1. Radial (Br) and vertical (-Bz) components of the magnetic field due to the disc currents. The magnitude of the dimensionless 
dipole field at the centre is 30 for models of type I and 60 for models of type II 
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Figure 2. Magnetic field lines for (a) models of type I and (b) models of type II. 
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Figure 3. Angular velocity f2 for ipo = 0.03 and /Bq = 0.1 (solid line), Vo = 0.03 and /3o = 0.04 (dotted line) and ipo = 0.06 and /3o = 0.1 
(dashed line). 

where — Bz/Bq. We first fix the value of |sm|. It is reasonable to assume that in an accretion disc the magnetic support 
due to the internal field is larger in the inner parts of the disc than in the outer parts. This is indeed the case if the magnetic 
field in the disc is advected radially due to the accretion (Lubow et al. 1994, Reyez-Ruiz & Stepinski 1996, Agapitou & 
Papaloizou 1996). In addition, the dipole increases the magnetic support in the central region of the disc. For these reasons 
we take |sm| to be a decreasing function of x. E is then calculated from |sni| and adjusted in order to make sense physically, 
n is then deduced from the dimensionless form of the radial equilibrium (^) 

fi' = fiic-/3o^, (30) 

where has also been used. As mentioned above (section 2.1), the inner parts of the disc where the magnetic field lines 
are linked to the dipole corotate with the dipole. The dipole flux is negative whereas the disc internal flux is positive. In the 
inner parts of the disc the total flux is then negative and in the absence of singular points the field lines are connected to the 
dipole. In contrast the total fiux in the outer parts of the disc is positive and the field lines are not linked to the dipole. We 
then fix the angular velocity f2 to be constant in the part of the disc where the total flux is negative, equal to its value at the 
point where the fiux vanishes. Given this new Cl profile, we then recalculate E using (^o|). 

We have performed calculations for disc models with a magnetic field configuration of the type described above. For all 
the models the inner disc radius Xin = 0.1 and the outer disc boundary is at Xo = 1. Figures ^ to ^ show respectively the 
equilibrium profiles of Cl, E, Sm and used for models of type I and II with (3o = 0.1 or 0.04. 

From Fig. ^ we see that the inner parts of the disc are dominated by the dipole field. Indeed /3o, which controls the 
strength of the disc magnetic field, does not affect the magnetic support close to the inner edge of the disc. 



5 NORMAL MODE CALCULATIONS 

We have performed global normal mode calculations for the disc models described above. Both axisymmetric modes and 
modes with small values of m have been considered and unstable spectra found. 
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Figure 4. Surface density S for i/iq ■ 
(dashed line). 



: 0.03 and /3o = 0.1 (solid line), Vo = 0.03 and /3o = 0.04 (dotted line) and Vo = 0.06 and /3o = 0.1 



We solve equation (^9|) by dividing the radial interval [a;in,a:o] into a grid of Ur points at positions (xi)^^^ ^ with a spacing 
Axi = Xi+i — Xi. The disc boundaries are at xi — Xin = 0.1 and — Xo = 1- 

We approximate the integral in the expression for by a sum of functional values at the points (xi)^^-, „ . Since the 



integrand involves through B^"^ (see (Q)), the discretised form of equation ( ^ ) is a system of equations of the form 



(a + mf2i) — 



i = 1 . 



(31) 



where the subscript i (or j) denotes the value of the function at the point Xi (or Xj) and A is a matrix which depends on the 
characteristics of the disc. If m = 0, ( pl[ ) is a x rir eigensystem with eigenvalue a and eigenfunctions £,z,i,i ~ 1 . . . n^. For 
m non-zero, we set 



2mQ, 



The system ( pl|) is then equivalent to a 2nr x 2nr eigensystem with eigenvalue a and eigenfunctions = 1 . . . and 

Ui,i — 1 . . .Ur - The eigenvalues are calculated numerically using the QR algorithm for real Hessenberg matrices given by Press 
et al. (1986). Once a is found, the solution of the system ( |3l| ) gives ^z,i,i = 1 . . . n^. 

Table ^ summarises the characteristics of the disc models considered. For a mode to be unstable we require that lm{a) < 0. 
The table also gives for each model the eigenvalues a for the unstable modes found. The real part of a is the frequency of the 
mode and the imaginary part relates to its growth rate. In the column which contains rir, we have indicated whether the grid 
is uniform (u) or non- uniform (n). 

The grids corresponding to rir = 864, 1175 and 1301 are non- uniform. They are characterised by a step Axi between a; = 0.1 
and 0.3, Ax2 between x = 0.3 and 0.8, and Axs between x = 0.8 and 1. The values of Axi, Ax2 and Axs are given in table ^ 
where Axo = 0.9/399. 

We performed a test for model I2a by setting i/)ext = 0. In this case with no external field, we confirmed that the disc has a 
zero frequency rigid tilt mode {^^ oc r) which is the mode having the lowest frequency. When an external field is added, the 
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Figure 5. Magnetic support Sm for ipo 
/3o = 0.1 (dashed line). 



0.03 and fio = 0.1 (solid line), Vo = 0.03 and fio = 0.04 (dotted line) and tpo = 0.06 and 



Table 1. Characteristics of the disc models and values of cr for unstable modes. 



Model 



cr (the unresolved modes are followed by *) 



I4b 



15 
II 



0.03 0.1 799 (u) (0,-37.44) (0,-17.31) (0,-7.62) 

— — 1 799 (u) (-8.77,-37.22) (-8.77,-17.80) (-8.77,-7.31) 

— — — 1175 (n) (-8.79,-36.90) (-8.79,-17.47) (-8.79,-6.72) 

— — 2 799 (u) (-17.54,-37.16) (-17.54,-18.21) (-17.54,-7.41) 

(-5.07,-0.0013)* (-4.49,-0.0015)* (-4.02,-0.0032)* 

— — — 864 (n) (-17.54,-37.14) (-17.54,-18.20) (-17.54,-7.40) 

(-6.78,-0.0008)* 

864 (n) (-26.31,-37.13) (-26.31,-18.60) (-26.31,-7.70) 

(-17.16,-0.0029)* (-7.70,-0.0019)* (-5.83,-0.0035)* (-5.39, 

— — — 1301 (n) (-26.31,-37.13) (-26.31,-18.60) (-26.31,-7.70) 

(-11.26,-0.0013)* (-9.17,-0.0012)* (-6.53,-0.0026)* (-6.30,- 
(-5.83,-0.0024)* (-5.40,-0.0037) 

— 0.04 1 799 (u) (-9.00,-37.07) (-9.00,-17.60) (-9.00,-6.88) 
0.06 0.1 1 799 (u) (-5.93,-40.33) (-5.93,-28.16) (-5.93,-23.06) 

(-5.93,-17.41) (-5.93,-12.95) (-5.93,-6.91) 



II 
I2a 
I2b 
I3a 

I3b 



I4a — — 



-0.0019) 
0.0022)* 



Table 2. Characteristics of the non uniform grids 



rir 


Axi 


Ax2 


Ax3 


864 


Axo/2 


Axo 


Axo/ 4 


1175 


Azo/8 


Axo/2 


Axo/2 


1301 


Axo/2 


Axo 


Axo/S 
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Figure 6. Square of the epicyclic frequency, , for i/iq 
ipo = 0.06 and /3o = 0.1 (daslied line). 



0.03 and /3o = 0.1 (solid line), Vo = 0.03 and f3o = 0.04 (dotted line) and 



equivalent mode is a modified rigid tilt mode, as shown on Fig. ^. Since the dipole field is significant only in the inner parts 
of the disc, this mode takes on the character of a rigid tilt mode in the outer parts. 

For m = and m = 1, all the unstable modes are well resolved. For higher values of m, in addition to well resolved unstable 
modes we get some poorly resolved weakly unstable ones. The number, frequency and growth rate of these modes depend on 
the grid resolution. However, in all cases they have growth rates several orders of magnitude smaller than those of the well 
resolved modes. Thus, although their reality is questionable, they are not important, and from now on we shall consider only 
the spectrum consisting of the well resolved unstable modes. 

These modes are confined in the inner parts of the disc where < (see Fig. ^ . This is in agreement with the local dispersion 
relation which predicts that instability ensues when the condition ( |l^ is satisfied. The frequency of these modes is —mClc, 
where flc is the angular velocity in the inner parts of the disc, and their growth rate depends only weakly on m. For each m, 
the most unstable modes have a growth rate significantly larger than their frequency indicating dynamical instability. The 
number of modes in the unstable spectrum increases with the strength of the magnetic support in the inner parts of the disc, 
as does the growth rate of the most unstable mode. 

For ^0 = 0.03 there are 3 unstable modes, with respectively 0, 1 and 2 nodes in the real part of ^2, whereas for ^po = 0.06 
there are 6 modes, with the number of nodes varying between and 5 (the smaller the growth rate, the larger the number 
of nodes) . Fig. ^ and Fig. ^ show the real part of ^2 in the inner parts of the disc for models I2a and II respectively. In all 
cases the imaginary part of is very small compared to its real part. We note that the characteristics of these modes do not 
depend on the resolution. 

Varying Po hardly changes the characteristics of the unstable modes. As mentioned above, the inner parts of the disc, where 
the modes are confined, are indeed dominated by the dipole field where conditions are insensitive to j3o. 

In the equilibrium models to which the above calculations correspond, the internal vertical field vanishes at some location 
in the outer parts of the disc. To check whether the results depend on the unlikely presence of this O-point, we have re-run 
model I2a and model 13 with rir — 1301 and the outer boundary at 2; = 0.7. The results are very similar to those obtained 
above, the minor differences coming from the fact that the rotation and the density profiles are a bit different in the two 
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Figure 7. Re(5z) for rpo = 0.03, (3o = 0.1, m = 1 and rir = 799 (model I2a). The mode represented is the modified rigid tilt mode and 
has a = (-0.025,0). 



cases. This is consistent with the fact that the internally confined unstable modes are almost independent of the structure of 
the outer parts of the disc. 



6 DISCUSSION 

In this paper, we have studied the stability of a magnetized accretion disc to disturbances perpendicular to its plane (bending 
modes). At equilibrium, the disc is permeated by both an internally produced poloidal magnetic field and an external dipole 
field. The former arises from a toroidal current in the disc, the latter is supposed to originate from a magnetized central 
star. In the important inner regions of the disc, the field lines are in a state of isorotation. We suppose that the uniformly 
rotating inner disc could be produced by material diffusing inwards through the action of interchange instabilities (Spruit & 
Taam 1990) when the mass injection rate into the disc is high enough to enable that to occur with ultimate accretion onto and 
spin up of the central star. We have neglected a possible toroidal component of the magnetic field which would be produced 
through differential rotation if poloidal field lines connected points with different angular velocity. Since the disc is supposed 
to be infinitesimally thin, radial pressure gradients have been neglected. 

A local stability analysis leads to the m independent condition for instability < (see (^)), where Km is the modified 
epicyclic frequency defined by the relation (]2^). In the inner regions of the disc, where the dipole field predominates over the 
internal one, this condition is satisfied if the magnetic field provides enough support against gravity (see ([l^)). We note that 
even though uniform rotation acts in favour of the instability, bending modes in an entirely differentially rotating disc may 
also be unstable. An external dipole is not absolutely necessary for instability to occur. 

For axisymmetric modes, the existence of a variational principle leads to the rigorous global criterion ( po| ) which is equivalent 
to ( [l6| ) when local displacements are considered. Even though the razor-thin disc approximation has been used so far, the 
instabilities are not an ab initio artefact of this assumption. The criterion ( po| ) is the thin disc limit of the general variational 
principle of Papaloizou & Szuszkiewicz (1992) for axisymmetric modes. 

We have solved numerically the normal mode equation ( |l^ ) for specific disc equilibrium models. For low values of m we get a 
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Figure 8. Re(5z) for 4>o = 0.03, /3o = 0.1, m = 1 and rir = 799 (model I2a). The modes represented have Re((7) = —8.77 and 
lm{a) = —37.22 (soUd hne), Im((T) = —17.80 (dotted hne) and Im((T) = —7.31 (dashed line). Only the inner parts of the disc are 
represented. 



spectrum of well resolved dynamically unstable modes, which are confined in the inner parts of the disc where < 0. For a 
given m, the number A'' of these modes increases with the magnetic support. For these the tendency is that the number of nodes 
in the real part of the vertical displacement increases from to — 1. This is similar to what happens in a non-magnetized, 
self-gravitating, uniformly rotating disc (we have already mentioned in section 3.2 the similarity between the two situations). 
In that case. Hunter & Toomre (1969) have indeed shown that the normal modes of free oscillation are polynomials with 
increasing numbers of nodes. However, in the self-gravitating disc, the bending modes are always stable. 

We have found numerically that the pattern speed associated with the unstable modes is ~ fie, with Sic being the constant 
angular velocity in the inner parts of the disc. 

We have argued that the instability of axisymmetric modes can be thought of as an unstable interaction between the disc 
current and the central dipole. Instability arises from the fact that the energy of the dipole in the magnetic field produced by 
the disc currents decreases due to the perturbation. 

The determination of the outcome of these instabilities awaits a non-linear analysis. However, we comment that if such 
instabilities occur, a configuration where the inner regions of the disc are displaced from the equatorial plane of the central 
star may be possible. Then the accretion of the disc material along the dipole field lines will be facilitated and may occur 
preferentially onto one stellar hemisphere depending on the mixture of normal modes present. Non-axisymmetric instabilities 
may result in a quasi periodic light variation even when the disc angular momentum vector at large distances and the central 
dipole axis are aligned. The asymmetric magnetospheric accretion would result in an observational signature of the star-disc 
configuration. The periodic light variation of Classical T Tauri stars is mostly interpreted as rotational modulation of the stellar 
fiux by hot spots due to magnetospheric accretion (see, for example, Bouvier et al. 1995). If accretion occurs preferentially 
to one hemisphere, then depending on the orientation of the observer the periodic light variation will be observed. The non- 
axisymmetric bending of the disc plane would lead to hot spots being created on the stellar surface and to the rotational 
modulation of the stellar output, even when the axis of the stellar dipole is aligned with the axis of rotation of the accretion 
disc. 
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a; 




Figure 9. Re((z) for Vo = 0.06, /3o 
Im(a-) = -40.33 (solid line), Im(CT) = 
represented. 



= 0.1, m = 1 and rir = 799 (model II). The modes represented have Re((T) = —5.93 and 
-28.16 (dotted line) and Im((7) = —23.06 (dashed line). Only the inner parts of the disc are 
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Im((T) = —6.91 (dashed line). 



Hunter C, Toomre A., 1969, ApJ, 155, 747 

Jackson J.D., 1975, Classical Electrodynamics. Wiley, New York 
Kenyon S.J., Yi I., Hartmann L., 1996, ApJ, 462,439 
Konigl A., 1989, ApJ, 342, 208 
Konigl A., 1991, ApJ, 370, L39 

Konigl A., 1993, in Levy E.H., Lunine J.I., eds., Protostars and Planets III. Univ. Arizona Press, Tuscon, p. 641 
Lehmann T., Reipurth B., Brandner W., 1995, A&A, 300, L9 
Lepeltier T., Aly J. J., 1996, A&A, 306, 645 

Lubow S.H., Papaloizou J.C.B., Pringle J.E., 1994, MNRAS, 267, 235 

Montmerle T., Andr'e P., Casanova S., Feigelson E.D., 1993, in Lynden-Bell D., ed., NATO Advanced Worksop on Cosmical Magnetism. 
Kluwer, Dordrecht 

Montmerle T., Feigelson E.D., Bouvier J., Andr'e P., 1994, in Levy E.H., Luninc J. I., eds., Protostars and Planets III. Univ. Arizona 

Press, Tuscon, p. 689 
Ogilvie G., 1997, MNRAS, 288, 63 
Paatz G., Camenzind M., 1996, A&A, 308, 77 
Paczynski B., 1991, ApJ, 370, 597 

Papaloizou J.C.B., Lin D.N.C., 1995, ARA&A, 33, 505 

Papaloizou J.C.B., Szuszkiewicz E., 1992, Geophys. Astrophys. Fluid Dyn., 66, 223 
Papaloizou J.C.B., Terquem C., 1995, MNRAS, 227,553 
Popham R., Narayan R., 1991, ApJ, 370, 604 

Press W.H., Flannery B.P., Teukolsky S.A., Vetterling W.T., 1986, Numerical Recipes: The Art of Scientific Computing, Cambridge 

University Press 
Reyez-Ruiz M., Stepinski T.F., 1996, ApJ, 459, 653 

Shu F.H., 1984, in Greenberg R., Brahic A., eds.. Planetary Rings. Univ. Arizona Press, Tuscon, p. 513 
Spruit H.C., Taam R.E., 1990, A&A, 229, 475 
Spruit H.C., Taam R.E., 1993, ApJ, 402, 593 

Spruit H.C., Stehle R., Papaloizou J.C.B., 1995, MNRAS, 275, 1223 
Tagger M., Henriksen R.N., Sygnet J.F., Pellat R., 1990, ApJ, 353, 654 
Tayler R. J., 1987, MNRAS, 227, 553 



© 1997 RAS, MNRAS 000, [[-|l| 



18 V. Agapitou, J. C. B. 

Wu F., 1987, ApJ, 320, 418 
Yi I., 1994, ApJ, 428, 760 



Papaloizou and C. Terquem 



© 1997 RAS, MNRAS 000, 



